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We present a non-equilibrium quasiclassical formalism 
suitable for studying linear response ac properties of Joseph- 
son junctions. The non-equilibrium self-consistency equations 
are satisfied, to very good accuracy, already in zeroth itera- 
tion. We use the formalism to study ac Josephson effect in a 
ballistic superconducting point contact. The real and imag- 
inary parts of the ac linear conductance are calculated both 
analytically (at low frequencies) and numerically (at arbitrary 
frequency). They show strong temperature, frequency, and 
phase dependence. Many anomalous properties appear near 
4> = IT. We ascribe them to the presence of zero energy bound 
states. 



I. INTRODUCTION 

Quasiclassical theories are proven to be powerful for 
studying superconducting systems. They have been used 
by many researchers to study the properties of supercon- 
ductors in the Meissner [1] or vortex states [2-4] , as well 
as for surfaces [5,6] , point contacts [7-10] , or grain bound- 
aries between different superconductors [11-14]. The 
equilibrium quasiclassical theory was developed by Eilen- 
berger [15] and Larkin and Ovchinnikov [16] by integrat- 
ing out irrelevant small scale degrees of freedom from 
the Nambu-Gorkov Green's function formulation of BCS 
superonductivity [17]. It was later generalized to non- 
equilibrium by Eliashberg [18] and Larkin and Ovchin- 
nikov [19]. 

A major development in the numerical calculations 
in equilibrium quasiclassical theory was made after the 
introduction of Riccati-transformation by Schopohl and 
Maki [2,20]. The transformation changes the Eilenbcrger 
equations [15] into a set of decoupled non-linear dif- 
ferential equations, that can be integrated easily. In 
non-equilibrium, the presence of convolution integrals in 
the equations of motion for the Green's functions makes 
the formalism nontrivial. Nevertheless, a generalized 
version of the Schopohl-Maki transformations for non- 
equilibrium systems has been suggested by Eschrig et al. 
[21,22]. 

Josephson junctions arc important devices, not only 
due to their rich physical properties, but also for many 
applications, including sensitive magnetometers [23,24], 
ultrafast switching devices [25], qubit prototypes [26,27], 
etc. dc and ac properties of them have been the subject 
of extensive research [28-33] . Some of the investigations 
are based on the tunneling Hamiltonian approach [34], 
which provides a good approximation when the trans- 
parency of the junction is small (e.g. tunnel junctions). 
At large transparencies, which is the case for supercon- 



ducting point contacts or grain boundary junctions, mul- 
tiple Andreev reflections (MAR) take place [35]. The 
MAR theories work well when the biasing voltage is large. 
At small biasing voltages, the number of Andreev reflec- 
tions grows (~ A/ey, with A being the superconducting 
order parameter). Nevertheless, the formalism was ap- 
plied to the case of single channel superconducting quan- 
tum point contact with small biasing voltage [31]. Alter- 
natively, non-pcrturbative Hamiltonian method [10] and 
non-equilibrium Green's function method [33] were de- 
veloped. However, most of these theories are suitable 
only when the applied voltage is constant. An exception 
is Ref. [32], which studies a superconducting quantum 
point contact in the presence of an ac voltage, but only 
in adiabatic regime. Thus a theory capable of studying 
high transparency Josephson junctions with an ac biasing 
voltage at arbitrary frequency is still lacking. 

For equilibrium systems, the quasiclassical Green's 
functions theory has provided a convenient tool to study 
Josephson structures with arbitrary transparency and 
roughness of the junctions (see Ref. [13] and references 
therein) . A generalization of the formalism to the non- 
equilibrium case, may also provide a powerful and conve- 
nient method to study ac properties of such structures. 
In this article, we rewrite the theory of Refs. [21,22] in a 
form suitable for studying properties of a general Joseph- 
son junction in the presence of an ac voltage at arbitrary 
frequency. To our knowledge, no such theory exists. We 
apply the theory to the case of a ballistic point contact 
between two conventional (s-wave) superconductors. At 
low frequencies, we find closed analytical expressions for 
the real and imaginary parts of the ac conductivity. They 
agree very well with the numerical results, except where 
the low frequency expansion fails. Both quantities show 
strong temperature, frequency, and phase dependence. 
We observe anomalous behavior when the phase differ- 
ence across the point contact approaches tt. We relate 
that to the presence of the zero energy bound states. 

Section II introduces the formalism and formulates it in 
a form appropriate for studying Josephson systems. Sec- 
tion III is devoted to calculation of ac current through 
a ballistic point contact between two s-wave supercon- 
ductors. The low frequency analytical results are given 
in III-A. The numerical results, as well as a comparison 
with the analytical results, are presented in III-B. Section 
IV summarizes the main results. A detailed description 
of the theory and notations is given in two appendices. 
Understanding the appendices in great detail is not nec- 
essary for understanding the main body of the article and 
for application of the theory to other problems. 
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II. THE FORMALISM 

A convenient way to study non-equilibrium systems is 

to use Kcldysh Green's functions [36]. The quasiclassical 
approximation of the Keldysh formalism introduces Re- 
tarded, Advanced, and Keldysh Green's functions: g", 
a = R,A, K. The first two describe spectral distribution 
of the states of the system, while the latter has infor- 
mation about population of those states. They are all 
2x2 matrices, and functions of the Fermi velocity vp, 
quasiparticle energy e, position r, and time t. The exact 
definitions of 5" and their equations of motion arc given 
in appendix A [see Eqs. (A2) and (A6)]. The method 
we present here is a linear response treatment of g^. We 
consider a clean superconducting system in the absence 
of an external magnetic field. The coupling to the elec- 
tromagnetic field is via the vector potential A and scalar 
potential $. As will become soon clear, working with 
a gauge in which A = simplifies the calculations sig- 
nificantly. In such a gauge, $ is the only perturbation 
applied to the system, which we take to be small. 

Let us first consider the case of a uniform (bulk) su- 
perconductor. We introduce a gauge transformation (see 
appendix A for details) A i-^ A = e~*'''^A. Under such a 
transformation e5> e$ — (l/2)9t(50. (Throughout this 
article we use 7i = fcs = 1.) We choose 5(j) in such a 
way to exactly eliminate <& from the gauge transformed 
dynamical equations [Eq. (A6)]. Thus, 5(1) should satisfy 

f-* (1) 

which is the well-known Josephson relation [37]. The 
vector potential A is still zero after this gauge transfor- 
mation, because W54> = within the bulk. The gauge 
transformed equation of motion is therefore exactly the 
same as the equilibrium equation. As a result, A = Aq, 
or A = Aqc*'''^, where Aq is the equilibrium order pa- 
rameter. This means that the magnitude of the order 
parameter, [A], is independent of whereas its phase 
varies with i> via the Josephson relation (1). In other 
words, if the only perturbation to the system is through 
the time varying potential its only effect is to change 
the phase of the order parameter [38]. 

The above simple observation has very important con- 
sequence in our linear response formalism. One can write 
the original (gauge dependent) order parameter as 

A = Aoe'*^ « Ao(l + i^</.) = Ao + ^A, (2) 

where (5A is the non-equilibrium linear response correc- 
tion to the order parameter. In Fourier spacx;, the .loseph- 
son relation (1) becomes 54> = i2e$(w)/a;, which yields 

(5A(a;) = -— e$(a;), (3) 

This value indeed satisfies the self-consistency equation 
for the homogeneous superconductor [cf. Eq. (23)]. To 



ensure small perturbation requirement, one needs e<E> ^ 
Lo. Charge neutrality condition in the bulk is also satisfied 
automatically by using (3). This can be seen simply from 
the above gauge transformation argument, taking into 
account the neutrality of the unperturbed (equilibrium) 
system. 

We now consider a Josephson junction with an equilib- 
rium phase difference (j) and an ac voltage V = Vo cos tot 
across the junction. Far away from the junction, we take 
the phase of the order parameter and the scalar poten- 
tial on the left (L) and right (R) sides to be (t^B.,!. = 
and $R.L = ±(Vo/2) costjt, respectively. In general, the 
phase of the order parameter is space dependent. There- 
fore, performing the above gauge transformation will pro- 
duce a vector potential A = (c/2e)V(5(/> (cx jac/jc.buik, 
where jac is the ac current density in the banks and 
Jc,buik is the bulk critical current density), which inval- 
idates our arguments. However, in most practical sys- 
tems, jac jc ^ .7c, bulk, whcrc IS the Josephson cur- 
rent density. The corrections to Eq. (3) are therefore 
small and 0(jc/ic,buik) [39]. Thus one can still use (3), 
as a very good approximation, even when Aq is space 
dependent. This removes the necessity for an iterative 
procedure (the main obstacle in these types of calcula- 
tions) in order to self-consistently calculate SA, and sat- 
isfy charge neutrality condition (within the bulk) [21,22]. 
The equilibrium order parameter Aq, however, should 
be calculated self-consistently using the common itera- 
tive methods; convergence of such calculations is proven 
to be very good, especially when using the Matsubara 
technique [13]. 

A. Equilibrium Solution 

In equilibrium, the retarded and advanced Green's 

fimctions can be written in terms of the Riccati ampli- 
tudes ttg and 60 in a way very similar to the conventional 
method for the Matsubara Green's functions [13]. 

where 

The subscript "0" denotes equilibrium quantities. The 
Riccati amplitudes satisfy the following Riccati equations 
[41] 

vf ■ Va^ = 2ie"a^ - {a^ fA* + Aq, 
-vp.Vb^ = 2ie''b^-{b^fAo + A*, (6) 

where e" = e + is^r], with e and r] being the real and 
imaginary parts of the quasiparticle energy respectively. 
77 is the quasiparticle damping, related to the inelastic 
lifetime r of the quasiparticles by = 1/r [40]. Notice 
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that the scalar potential $ does not appear in the equi- 
librium equations. The boundary conditions are the bulk 
solutions of Eq. (6): 



an = 



Ao 



Oq — 



A* 

^0 



(7) 



where 0"=i/|Ao|^ — (e")^. The differential equations 
(6) and the boundary conditions (7) can be obtained from 
their counterparts in Matsubara formalism (see Ref. [13] 
for example), by changing a;„ — »■ — ie", where a;„ is the 
Matsubara frequency. 

To calculate the Riccati amplitudes at other points, 
one should define quasiclassical trajectories as straight 
lines in the direction of vf (see Fig. 1). and are 
obtained by integrating (6) in the direction of along 
the trajectory, starting from the boundary conditions (7) 
at —00. For and , integrations are taken in the 
opposite direction. The following symmetries exist for 
the equilibrium functions 



(8) 



It is therefore sufficient to calculate one of the sets of 
retarded or advanced functions. 

In equilibrium, the Keldysh Green's function is related 
to the retarded and advanced ones by 



(5? 



where 



= tanh ( I 
\2TJ 



(9) 



(10) 



takes into account the thermal distribution of the quasi- 

particles. Equilibrium current and charge densities are 
calculated by Eqs. (A17) and (A18), summing over all 
trajectories. 



approximation to the exact solution. All the differen- 
tial equations [Eqs. (B2) and (B4)] will then be linear in 
the time-varying parts, and Fourier transformation will 
be straightforward: no complicated convolution integrals 
arise. 

Let us introduce simplifying notation 



e±=e±-, (9o±=0o(e±)- 



(12) 



We define linear response Green's functions as 6g" = 
—Qq, where are the equilibrium Green's functions. 
Similarly, we introduce small corrections to the Riccati 
amplitudes (5a" = a" — ag and Sb" = 6" — 6q . The linear 
response Green's functions are then given in terms of the 
Riccati amplitudes by 



7" = -2s° 



5r = 2s° 



(l + a^+6;?+)(l+a^_6g_) ' 
5a« - 5b°'a'^^a'§_ 
(l+ag+&^+)(l + ag_6g_) ' 



(13) 
(14) 



for a = R,A. 

We also define an anomalous Green's function 5g^ by 

Sg"" = 5g''{T+-T-) + 5g"'T--5g^:F+. (15) 

Correspondingly, we introduce anomalous functions 5a^ 
and 5b^ [see Eqs. (B9) and (BIO)] which are related to 
the Green's functions through 



Sg 



(l + a«X+)(l + «o3^) ' 
Sa^a^_ + 6b^a^^ 



(16) 
(17) 



The differential equations describing 5a" and 5b" {a 
R, A, X) have general forms 



B. Linear Response Solution 

Generalization of the above Riccati-transformation to 
non-equilibrium is discussed in appendix B. The presence 
of the (8)-operations [see Eq. (A8) for definition] makes 
calculations nontrivial. However, significant simplifica- 
tion arises when one side of the ^-operation is an equi- 
librium quantity (and therefore time/frequency indepen- 
dent). More specifically, in frequency space we have 



Po(e) ® Q{e,u;) = Po [e + -) Q(e,c^), 

P(e,a;)®Qo(e) = P(e,w)Qo(e-|). (11) 

To take the advantage of this property, we use linear 
expansion. In other words, we assume that the pertur- 
bation to the system (i.e. $) is so small that the lin- 
ear expansion of the Green's functions provides a good 



-vf ■ VSb" = A^Sb^ + B", 
with ^'s and B's given by [41] 



(18) 



A" = 2ze-AUa^^+a^_), 

B" =5A + a^^a^_5A* - ie$(a[J+ - a^_), (19) 
I" = 2ze-Ao(6?+ + 6^), 

B" = -6A* - b^+b^_SA + ie<i>{b^+ - b'^_), (20) 

for the retarded (a = R) and advanced (a = A) func- 
tions, and 



A 



X 



ai,\A*+btAo, 



B^ = -a^^5A* + b^_5A - ie$(l + a^+b^_ 



(21) 
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A 



X 



b^^Ao+a^A*, 



b^ 



(22) 



for the anomalous ones (a — X). 5 A = A — Aq, is given 
by Eq. (3) and satisfies the self-consistency equation 

<5A(v^) = ^ de (VivF, V^)Sf^{V^))^,^ , (23) 

where V{vF,v'p) is the interaction potential, Np the 
density of states at the Fermi surface, and Cc the energy 
cutoff. The bulk boundary conditions for the amplitudes 
are 



5b°' 



(24) 



To assure stability of the integrations, Eqs. (18) should be 
integrated along the trajectory in the direction of vj? for 
(5a^, bb^ , and , but in the opposite direction for (56^, 
ba^ , and bb^ (see Fig. 1). Having the Green's functions, 
the non-equilibrium correction to the current density is 
given by 



rfe(vF Tr[f3(5?^]) 



and the charge density is found from 

bp = eNp \ -2e^ + i / '^'^ i^^i^d^]) 



(25) 



(26) 
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FIG. 1. Two superconducting regions connected via an ori- 
fice. The dashed line shows the quasiclassical trajectory. The 
arrows indicate the directions of integration. 

that the linear expansion provides a good approximation 
{eVo < |Ao|,L^). 

To find the current response, /, of the system, we calcu- 
late the current density bjz at the orifice, using Eq. (25), 
and then integrate it over the area S of the orifice. In 
Fourier space, I{lo) can be a complex number. The real 
part of it describes the dissipation of the system, while its 
imaginary part gives information about inductive or ca- 
pacitive behavior of the system. The linear admittance of 
the system is defined by F(u;) = /(tj)/Vb (= ^/Z, where 
Z is the impedance of the contact). 

We proceed with the calculation of the current in two 
different ways. First, we find analytical results in the 
regime of small w and rj. We then provide the results 
of full numerical calculation and compare them with the 
analytical ones. 



III. AC JOSEPHSON EFFECT IN A 
SUPERCONDUCTING POINT CONTACT 



A. Low Energy Analytical Results 



A ballistic superconducting point contact is probably 
the simplest system that the method described here can 
be applied to. It brings extra simplicity because even the 
equilibrium solutions can be found non-self-consistently 
[7]. Nevertheless, the system shows rich and nontrivial 
physical behavior. Let us consider an orifice between 
two conventional (s-wave) superconductors (Fig. 1), di- 
mension of which much smaller than the inelastic scat- 
tering length and coherence length of the superconduc- 
tors {d <^C Zt-=i^ft, ^o=i'F/7r|Ao|). We assume perfect 
transparency at the contact, although generalization to 
arbitrary transparency is straightforward [13,22]. We 
take the equilibrium order parameter to be constant 
on both sides of the contact, Al,r — lAole*"*"--^, with 
0L,R — i0/2, where </> is the phase difference between 
the two sides. Here, L (R) denotes the left (right) side of 
the contact. 

We also take i5A to be constant on each side of the con- 
tact, given by Eq. (3). The scalar potential $ is taken 
to be $L,R = ±(Vb/2) coswt on the left and right sides 
respectively, Vq being the amplitude of the potential dif- 
ference. It is taken to be a real number and so small 



In the case of a point contact, it is not difficult to ob- 
tain analytical results. Since the superconductor is ho- 
mogeneous everywhere except near the contact, the solu- 
tion to the functions a and b at the contact is almost equal 
to their bulk values. More specifically, a^(0) = a^(— oo), 
6^(0) = 6^(-foo), etc. From now on we drop the ar- 
guments and just write a^, 6^, etc., keeping in mind 
that what we mean is the values at the position of the 
contact. Substituting these values in the corresponding 
equations, one can obtain analytical expressions for the 
current. The exact expression is rather complicated and 
does not give any more insight than the numerical re- 
sults. It, however, can be significantly simplified in low 
energy regime. 

Here, we calculate the contact admittance in the 
regime 77,0; <C |Ao|,T (but of course w » eVb)- Let 
us first introduce the following parameterization 



e" = |Ao|cos7", a^R.A 



(27) 



where 7" is a complex number. We therefore find SI" — 
I Aq I sin 7". This choice of notation significantly simpli- 
fies the form of ao and 60 • For instance from (7), taking 
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the phase of the order parameter to be +(j>/2 (left side of 
the contact), one finds 



an = 



-le 



(28) 



where = 7^ — 0/2. Similarly, has to be calculated 
on the right side of the contact (where the phase is — (?i/2) 
and turns out to be exactly equal to aQ. In general, one 
can write 



(29) 



First notice that 1 + a^6^ = 1 — q-'^^^''^'i± vanishes as 
5^l'^ = Sj"{e'^) (or as 7^ 0/2)- Because of such 
expressions in the denominators of Eqs. (13) and (16), 
the most important contribution to those equations must 
come from points close to e = eo{(p) for which ^7" <C 1. 
Here 



eo(0) = |Ao| cos ■ 



(30) 



(±eo are indeed the energies of the Andreev bound states 
in the contact [36].) We can therefore focus only on those 
points. Expanding the numerators of (13) and (16) up 
to first order in 6j", (i.e. first order in lu and rj), we find 



, A el/ 



1 



eV 



1 

1 



+ 



+ i 



+ i 



+ 



2^0 



— z 



where 



Oo(<^) 



2f2o 

2O0 
1 



1^0 1 sin-. 







(<57f 


,57«). 










1 — ^ 




57f 





(31) 



(32) 



Let us write e = eo + e' (e' -C |Ao|). For 0^0 (we will 
discuss the <p = case later), one can write 



Therefore 
1 



p 



e' ±Lu/2 



— ins 



.(. 



± 



(33) 



(34) 



where P gives the principal value integral when integrat- 
ing over e'. Because of symmetry, the principle value 
integrals are negligible after integration. We therefore 
write 



(35) 



On the other hand 

1 



(57^(57^ 

Similarly 

1 



CO 



"0 



■uj/2 + is°'rj e' + uj/2 + is°"q 



<^7+57f w + 2iri \e' - uj/2 - irj e' + uj/2 + i-q^ 



LO + 2i?7 



Substituting (35)-(37) into (31) and using (15), one can 
calculate Sg^ and thereby the total current using (25). 
Taking the integral over e' in (25), expanding the result- 
ing hyperbolic tangents around eo (to the first order in 
w/T <C 1) and keeping only the leading order terms, we 
find the admittance to be 



y(a;) = 



■N 



27] 



2T \L0 + 2ir] 



sech^ 

2T 



+ ieo tanh ^ 



(38) 



where Rn ~ 2/ e^vpNpS is the normal (Sharvin) resis- 
tance of the point contact. Equation (38) agrees with 
the result obtained by Averin and Bardas [32] in adia- 
batic regime [42] . The quasiparticle conductance is given 
by the real part of (38): 





Gn (w2 + Arf)T 



,2 ^0 
sech — , 
2T 



(39) 



where Gat = 1/Rn. Notice that the right hand side of 
Eq. (39) vanishes at = 0. This however is the point 
at which the linear expansion (33) fails. One therefore 
expects the terms neglected in (39) to dominate G{(l) = 
0). Similarly, at small T, Eq. (39) is exponentially small 
(except at (/) = tt). Since in such a regime, lo/T becomes 
large, the expansion of the hyperbolic tangents in powers 
oiu/T will be invalid. Therefore, deviation from (39) at 
low T is expected. In the next subsection, we will see 
such deviations in the numerical results. 

The imaginary part of the admittance is also impor- 
tant because it provides information about the inductive 
or capacitive behavior of the junction. Notice that at 
(f) = 0, the first term in (38) vanishes but the second 
term survives. Thus unlike the conductance, Im(F) does 
not vanish; it rather stays finite and behaves purely in- 
ductively (~ I/ojL). Similarly, near T = 0, the second 
term in (38) dominates, resulting in a finite (again induc- 
tive) Im(y ) . One therefore expects that the leading order 
expansion provides good approximation. The exception 
is at ^ = TT where eo = and thus the second term in 
(38) vanishes. The higher order terms therefore play im- 
portant role in such a case. In the next subsection, we 
observe this behavior by comparing with the numerical 
results. 
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B. Numerical Results 

We now present the results of our numerical calcula- 
tion. The value of |Ao| is calculated directly from the 
BCS gap equation [34]: 



1 = A 



de 



|Ao| V« 



= tanh ( — I 



(40) 



where the dimensionless coupling constant A is chosen in 
such a way to give Aq^O as T^T^- Near T = 0, one finds 
|Ao| « 1.75Tc. All the energy scales (T, e, w, eVb, 77, etc.) 
are normalized to Tc and y{iij) to Gat. In all calculations, 
we take -q = 0.01 and tc — 20. 



Numerical 
Analytical 



CO 



FIG. 2. Linear conductance G (normalized to Gjv) as a 
function of ui, for </> = 37r/4 at T = 0.1. (u and T as well as 
other energy scales are normalized to Tc.) 

Fig. 2 compares the result of numerical calculation of 
G ai (f) — in / A with the analytical result obtained from 
Eq. (39). As expected, the two curves overlap at low fre- 
quencies but deviate at larger ld. Around w = 0, there 
exists a sharp peak corresponding to the Lorentzian lj- 
dependence in (39). At larger frequencies, a second peak 
appears in the numerical curve which is absent in the ana- 
lytical one. The peak clearly results form the higher order 
contributions which were neglected in derivation of (38). 
Fig. 3 displays numerical G-oj curves at different phase 
differences. At = tt (Fig. 3d), the sharp peak at w = 
has the largest value. At smaller phase differences, this 
peak becomes less pronounced and eventually disappears 
at = 0, as it should according to (39); the equation 
actually predicts zero conductance, therefore the small 
conductance in Fig. 3a is completely due to the terms 
neglected in (39). The second peak, however, appears at 
uo = |Ao| w 1.75 for cj) = Q (Fig. 3a), and moves towards 
w = as — > TT. It is easy to see that the peak is always 
at a; = eo(0)i i-e., the energy of the Andreev levels. 

Fig. 4 displays the conductance G as a function of the 
phase difference across the contact for different frequen- 
cies. It is clear from the figure that G is strongly phase- 
dependent. Especially, it is sharply peaked close to 4> = n 
[for bj = 0.01, it is more than five orders of magnitude 
larger than G(<^=0), see Fig. 4b]. The strong conduc- 
tance can be attributed to the existence of zero energy 




FIG. 3. Linear conductance as a function of lj at T = 0.1. 
The phase differences are (a) = 0, (b) = 7r/2, (c) 
(j) = 37r/4, and (d) (f> = n. 
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FIG. 4. (a) Linear conductance G as a function of the phase 
difference at T = 0.1, for different frequencies u. (b) The 
same data on logarithmic scale. The solid lines correspond to 
the analytical results. Here (and in the figures that follow), 
the legend is common between (a) and (b). 



Andreev bound states (ZBS) [i.e. €o{4'=tt) = 0]; they 
provide large density of states at zero energy. A compar- 
ison with the analytical results is shown, in logarithmic 
scale, in Fig. 4b. As expected, the agreement between 
the two calculations at w = 0.01, is good near (j) = n but 
they deviate as — > 0. The curves however overlap less 
at higher frequencies. Especially, at u; — 1 they show 
completely different phase dependence. 

The temperature dependence of the linear conductance 
is presented in Fig. 5a. All the curves join at G = 1 (or 
G = Gm before normalization) as T^Tc- This indeed 
is expected, because at T = Tc the superconductor be- 
comes normal. It is clear from the figure that the conduc- 
tance behaves completely differently at = tt compared 
to other phase differences. At = tt, the conductance 
grows with lowering the temperature and exhibits a 1/T 
dependence in agreement with (39). This type of l/T be- 
havior also exists in the dc Josephson current at = tt, 



6 



(a) 
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FIG. 5. (a) Linear conductance as a function of tempera- 
ture for different phase differences at u} = 0.1. (b) The same 
data plotted in logarithmic scale. Solid lines are analytical 
results. Notice that at <^ = 0, Eq. (39) gives G = 0, therefore 
no analytical curve is shown in the figure. 



(b) 




0.02 0.04 0.06 0.08 0.1 0.02 0.04 0.06 0.08 0.1 

FIG. 6. (a) Conductance G as a function of damping rate 
r] a.tuj = 0.005. (a) T = 0.01, (b) T = 0.5. 



and is associated with the current carried by the ZBS. 
For (j) = 0, 7r/2, and 37r/4, and at intermediate tempera- 



tures, the Hnear conductance behaves as G 



-eo(0)/T 



in agreement with Eq. (39). This form of suppression re- 



sembles the thermal activation behavior (G ■ 



-A/T 



) in 



tunnel junctions [34]. At lower temperatures however, a 
deviation from such a behavior occurs. To examine this 
more carefully, and also to compare with the analytical 
results, we have plotted the same graph in logarithmic 
scale in Fig. 5b, adding to it the analytical curves (solid 
lines). Except for the = case [where (39) vanishes], 
the agreement between the numerical and analytical re- 
sults at intermediate temperatures is very good. At low 
T, on the other hand, the numerical curves show satura- 
tion. The crossover temperature to the saturation regime 
is proportional to eo(0) and is almost independent of uj 
(not shown in the figure). 

Such a saturation does not occur in the analytical 
curves (naturally, was also not predicted in Ref. [32]), 
and is clearly a higher order property. As we mentioned 
before, the leading order result of (39) vanishes as T ^ 0, 
therefore the only remaining contribution will be the 
higher order terms. To see this explicitly, in Fig. 6 we 
have plotted G versus at a low frequency (w = 0.005) 
and for two different temperatures. One immediately no- 
tices significant difference in the ?7-dependence between 



the two cases. At high T (Fig. 6b), the conductance 
follows 1/77 dependence in agreement with (39). At 
T = 0.01 (Fig. 6a), on the other hand, all three curves 
show linear dependence on rj, which is obviously higher 
order than l/rj. Physically, the residual conductance is a 
result of the overlap of the midgap states, broadened by 
finite rj, at zero energy. Increasing 77, increases the density 
of states at zero energy and therefore the conductance. 
In reality, 77 is also temperature dependent and in s-wave 
superconductors, it vanishes at T = 0, and so does G. 




FIG. 7. (a) The imaginary part of the admittance as a 
function of at T = 0.1. (b) The same data multiplied by ui 
together with the analytical (solid) curves. 

Fig. 7a shows the frequency dependence of Im(y) for 
different phase differences. Except for the (j) = tt curve, 
the other curves seem to show a 1/w form (inductive be- 
havior). To see this more clearly, we have plotted ujlm{Y) 
as a function of uj in Fig. 7b. The first three curves are 
almost constant, confirming the 1/uj dependence. The 
curve at (/) = TT, on the other hand, exhibits completely 
different behavior. In Fig. 7b, we also present the ana- 
lytical curves corresponding to Eq. (38). The agreement 
between the numerical and analytical curves is good at 
low frequencies. At higher frequencies, all curves deviate 
from their analytical counterparts, as expected. Espe- 
cially, for the case oi (p — tt the discrepancy between the 
two curves is significant. 



E 
3 







-6 
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(i) = 0.01 

(0 = 0.1 

(0 = 2 



\ 



0.2 



0.4 0.6 
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FIG. 8. u) lm{Y) as a function of phase difference for dif- 
ferent frequencies at T = 0.1. The thin lines correspond to 
the low- frequency analytical results. 

To understand this better, we have plotted wlm(F) 
versus the phase difference (/) for different frequencies in 
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Fig. 8. The thin curves are plotted using Eq. (38). For 
uj = 0.01 and 0.1, the curves overlap and agree quite well 
with the analytical ones over a wide range of cf). Near = 
TT, the curves separate and deviation from the analytical 
results become more evident. At high frequency [uj — 2) 
on the other hand, the deviation already exists at </) = 
and increases as — > tt. 




FIG. 9. Im(y) as a function of temperature for different 
phase differences at uj = 0.1. The solid curves show the cor- 
responding analytical results. 

The temperature dependence of Im(y), for different 
values of the phase difference, is plotted in Fig. 9. All 
the curves meet at lm(y)=0 as T ^ Tc. The agreement 
with analytical results is very good at high T. At lower 
temperatures, the (j) = tt curve behaves completely dif- 
ferently than the other curves and deviates significantly 
from the analytical curve. This, as we mentioned be- 
fore, is a result of the breakdown of the small frequency 
expansion at low temperatures. 



IV. CONCLUSIONS 

We have presented a microscopic formalism for cal- 
culating ac properties of Josephson junctions. The 
method is based on linear response treatment of the 
non-equilibrium quasiclassical Green's function theory of 
superconductivity and uses a generalized form of the 
Riccati-transformations. The self-consistency equation 
for the linear response part of the order parameter, as 
well as the charge neutrality condition within the bulk, 
is satisfied to a very good approximation, with no need 
for numerical iteration. 

We successfully applied the method to the case of a bal- 
listic superconducting point contact and obtained non- 
trivial results for linear conductivity of the junction both 
analytically and numerically. We noticed strong temper- 
ature, frequency, and phase dependence in the real and 
imaginary parts of the ac conductivity. In particular, 
we found the conductance to be many orders of magni- 
tude larger at (j) — n than at smaller phase differences. 
This is a result of the influence of the zero energy bound 
states on the quasiparticle conductance. The agreement 



between the analytical and numerical results is very good 
at low frequencies. The exceptions happen near (f> = 
for G, and (j) = n for Im(y), where the leading order con- 
tributions vanish or become comparable to the neglected 
terms. The discrepancy becomes more pronounced at 
low T or high uj, where the validity of the leading order 
approximation become questionable. 

Experimentally, superconducting point contacts have 
been realized using techniques such as: scanning tun- 
neling microscopy [44], mechanically controllable break 
junctions [45,46], superconductor-2 dimensional electron 
gas-superconductor junctions [47,48], etc. Subgap struc- 
tures were observed [49] and nice measurements of trans- 
mission coefficients of individual quantum channels were 
performed by Scheer et al. [50]. Unfortunately, phase- 
dependent measurement of the conductance is difficult 
and the only available data (to our knowledge) is those 
by Rifkin and Deaver [51]. They found a strongly phase 
dependent conductance, in qualitative agreement with 
our results. More experimental research is necessary to 
confirm the predictions of the present work. 

In this article, we have only considered a contact with 
perfect transparency. The method, however, is general 
and applicable for arbitrary transparency. One only 
needs to take into account appropriate (e.g. Zaitsev) 
boundary conditions at the contact. A solution to the 
Zaitsev boundary conditions, suitable for the present 
calculation method, is given in Ref. [22]. The method 
proposed in this article is also applicable to other sys- 
tems such as grain boundary junctions between uncon- 
ventional superconductors, which is the subject of a sep- 
arate publication [52]. 
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APPENDIX A: QUASICLASSICAL KELDYSH 
GREEN'S FUNCTIONS 

The quasiclassical Green's function [43] g{vF,e;r,t) is 
a 2 X 2 matrix, every element of which also a 2 x 2 matrix, 
and a function of the Fermi velocity vp, quasiparticle 
energy e, position r, and time t. We represent g as 







(Al) 



where the matrices g^, g^, are the quasiclassical 
retarded, advanced, and Keldysh Green's functions in 
Nambu-Gorkov representation, respectively: 
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cR,A 



9 — \ j:R,A^ _gR,A'\ 



(A2) 



The f-operation performs the fohowing transformation 

0\wF, e; r, t) = 0{-Vf, -e; r, t)* . (A3) 

In frequency domain it also changes lu to —lu. The re- 
tarded and advanced Green's functions carry information 
about the energy spectrum of the electronic states of the 
system, while the Kcldysh ones have information about 
occupation of those states. 

The following symmetries hold for the Green's func- 
tions: 



and 



g^ = -5^*, 



5^ = 5^*, 



f 



f 



K* 



(A4) 



(A5) 



In frequency space, these give g^{uj) = —g^{—iv)*, etc. 

The equation of motion that describes the time evolu- 
tion of g is written as 



i (^e - -vf ■ a ) -fs - a + ie$i, g 



accompanied by the normalization condition 
g^g = l, 



:0 (A6) 



(A7) 



where 3> and A are the scalar and vector potentials re- 
spectively, and [A, B]^ = Ai^B-Bi^A, with 

{A ® B){e, t) = ei^^^^^-^^^?)A{e, t)B{e, t). (A8) 

The product (A8) is associative and satisfies {A® B)'^ = 
At (g) St, but {A O B)* = (g) A*. We also have 



T3 



fa 
fs 



A 
A 



where the 2x2 Pauli matrix % and A are 

A: 



T3 



1 

-1 



A 
At 



(A9) 



(AlO) 



with A being the superconducting pairing potential. The 
constant e in (A6) is the absolute value of the elec- 
tronic charge, and c is the speed of light. In equilib- 
rium, Eq. (A6) reduces to the Eilenberger equation (see 
Ref. [13] for example), by changing e itOn, where a;„ 
are the Matsubara frequencies. 

A gauge transformation is defined by the following si- 
multaneous operations: 



2e 

A A - 2^Vx, 



(All) 
(A12) 
(A13) 

(A14) 



It therefore takes the phase of the order parameter (j) to 
(j) + X- 111 tlie discussion of section II, we have chosen 

X = -H- 

The Green's functions can be used to calculate physical 
quantities. The quasiparticle density of states is given by 



A^(e) = ^(Tr[f3(5«-5^)]) 



(A15) 



where Np is the density of state of electrons at the Fermi 
surface and (...) denotes averaging over wp- The pairing 
potential satisfies the following self-consistency relation: 



A(vf) 



Nf 



fie(y(vF,v^)/^'(v^)) , , (A16) 



where is the off-diagonal part of and V{wf,'v'f) 
is the interaction potential. Furthermore, the current 
density is given by 



J = 



eNp 



de{vF Tr[f35^]), 



(A17) 



and the charge density by 



p = eNF{ -2e$ + - / de {Tr[g'']) ] . (A18) 



APPENDIX B: GENERALIZED 
RICCATI-TRANSFORMATION 

Riccati-transformations [20] are proven to be very use- 
ful tools for numerical calculations of equilibrium prop- 
erties of superconducting systems. For non-equilibrium 
systems, however, the presence of the (^-operators makes 
the formalism nontrivial. Nevertheless, a generalization 
of the standard transformation to the non-equilibrium 
case is possible. It is common to define Riccati ampli- 
tudes a" and b", where a=R, A for the retarded and ad- 
vanced Green's functions respectively. They are related 
to the corresponding Green's functions by [21,22] 

5" = s"(l + a" 6")-i (1 - a" (g) 6"), 

/" = s"(l + a" 5")-^ (2a"), (Bl) 

where s" is defined in Eq. (5). Here we define the inverse 
operation by A = ^0 = 1. Eqs. (Bl) already 
resemble their counterparts in the standard Matsubara 
formalism [20] . It is straightforward to show that for A = 
these functions satisfy Riccati-type equations given by 
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ViT • Va« = 2iea" - a" A^' O a« + A + [ie<E>, a"]^ , 
-VP . V6" = 2ie6" - 6" (8> A (8> 6" + At - [ie$, 6"]^ (B2) 

In equilibrium, the ^-operation is replaced by a simple 
multiplication and these equations reduce to (6). 

It is also necessary to define other functions and 
[53], which are related only to the Keldysh Green's 
functions [21,22] 

= 2(1 + ® {a^ + a^(^b^(» h^) 

O (1 + O 6^)"S 
= 2(1 + a« ® (a-^ ® - a^" ® b^) 

® (1 + 6^ (g) a^)-\ (B3) 

and are governed by the following dynamical equations 
VF ■ Va^ = -dta^ - (g) (g) + (g) A (g) b^ 



[ie$. 



.^1 



-vf • Vb 



a 

I 

[ie^,b^]^. (B4) 



^ - - 6^ (g A (g) 6^ + 6^ (g) At (g) 



The functions a" and 6" are related, by the f-operation 
[Eq. (A3)], through 6"=a"t for R,A, K. In addition, 
the symmetries (A4) and (A5) require 



(B5) 



One can show that Eqs. (B1)-(B4) satisfy the dynamical 
equation (A6) together with the normalization condition 
(A7). 

In equilibrium, Eqs. (Bl) reduce to (4), and (B3) give 



50^ = 2 



K 



(l + ao«&o^)(l + aX)' 
- a^b^ 



(B6) 
(B7) 



(l + aX)(l + «X)' 

Satisfying (9), one finds 

= (1 + a«6^)jr, b^ = -{l+ 4b^)J^. (B8) 

Linear response equations (13) and (14) are obtained 
by expanding (Bl) to the linear order. To obtain (16) and 
(17), we expand (B3) to the linear order and introduce 



1 - ( -5g^\ 



(B9) 



through (15). We also define 5a^ and 5b^ in terms of 
the linear response 5a^ and 5b^ by 

5a^ = 5a^{J'+ - TJ) + 8aH^_T^ + bb^a^^T^, 

hb^ = 6b^{T- - T+) - 5b^a^_T^ ~ 5aH^+T+. (BIO) 

The differential equations (18) then follow directly from 
(B2) and (B4). It should finally be mentioned that the 
A's and B's in Eq. (18) are related by: A« = A°'\ B° = 
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